Consideremos la siguiente integral bidimensional: \[\begin{equation*} I = \int_{1}^{4} \int_{2}^{7} \frac{x^{2}y}{3} \,dx \,dy = \int_{0}^{1} \int_{0}^{1} 5 (5 u_{1} + 2)^{2} (3 u_{2} + 1) \,du_{1} \,du_{2} \end{equation*}\]
A continuación vamos a estimar el valor de \(I\) mediante el método de Montecarlo, tanto de forma directa como usando distintas técnicas de reducción de varianza. Para poder comparar más fácilmente los resultados obtenidos, los recopilaremos en una tabla al final de este documento.
Para cada método estimaremos también su coste en tiempo haciendo uso
de las herramientas proporcionadas por el paquete bench (en
particular, la función mark analiza el coste en tiempo y en
memoria de las expresiones proporcionadas, ejecutando cada una de ellas
un cierto número de iteraciones y devolviendo una tabla con distintas
medidas, entre ellas la mediana de los tiempos de ejecución de cada
iteración). De esta forma, podremos comparar la eficiencia de cada
método a la hora de estimar el valor de la integral. Para ello, es
fundamental usar la misma unidad de tiempo a la hora de estimar el
coste.
unidad_de_tiempo <- "s"
genera_vector_aleatorio <- function() {
runif(2)
}
g <- function(u) {
u_1 <- u[[1]]
u_2 <- u[[2]]
5 * (5 * u_1 + 2)^2 * (3 * u_2 + 1)
}
n <- 1e4
coste_directo <- bench::mark(
{
valores <- replicate(n, {
u <- genera_vector_aleatorio()
g(u)
})
},
iterations = 10,
time_unit = unidad_de_tiempo
)$median
estimacion_directo <- mean(valores)
varianza_directo <- var(valores) / n
eficiencia_directo <- 1 / (varianza_directo * coste_directo)
En primer lugar definimos los estratos y establecemos la forma de generar valores aleatorios dentro de cada estrato. Una manera simple de realizar un muestreo sistemático (es decir, de considerar estratos todos del mismo tamaño) es subdividir el intervalo \((0, 1)\) de forma independiente en cada dimensión.
genera_subintervalos <- function(numero_subintervalos) {
extremos_derechos <-
seq_len(numero_subintervalos) / numero_subintervalos
extremos_izquierdos <-
c(0, extremos_derechos[-numero_subintervalos])
data.frame(
min = extremos_izquierdos,
max = extremos_derechos
)
}
# Como en cada dimensión se va a considerar el mismo número de
# subintervalos, basta guardar la información una única vez
numero_subintervalos <- 5
subintervalos <- genera_subintervalos(numero_subintervalos)
# Los estratos vienen dados por todas las combinaciones posibles
estratos <- expand.grid(
seq_len(numero_subintervalos),
seq_len(numero_subintervalos)
)
cantidad_estratos <- nrow(estratos)
estratos$probabilidad <- 1 / cantidad_estratos
genera_vector_aleatorio_en_estrato <- function(numero_estrato) {
numero_estrato_u_1 <- estratos[numero_estrato, 1]
numero_estrato_u_2 <- estratos[numero_estrato, 2]
estrato_u_1 <- subintervalos[numero_estrato_u_1, ]
estrato_u_2 <- subintervalos[numero_estrato_u_2, ]
c(
runif(1, min = estrato_u_1$min, max = estrato_u_1$max),
runif(1, min = estrato_u_2$min, max = estrato_u_2$max)
)
}
Ahora replicamos el proceso de generar valores en cada estrato, en
una cantidad proporcional a su probabilidad, y aplicarles la función
g a cada uno de ellos.
n_estratos <- (n * estratos$probabilidad) |>
ceiling() |>
pmax(2)
coste_estratificado_proporcional <- bench::mark(
{
valores <- purrr::map2(
seq_len(cantidad_estratos),
n_estratos,
\(numero_estrato, n_estrato) {
replicate(n_estrato, {
u <- genera_vector_aleatorio_en_estrato(numero_estrato)
g(u)
})
}
)
},
iterations = 10,
time_unit = unidad_de_tiempo
)$median
Finalmente, estimamos el valor de la integral, la varianza de esa estimación y la eficiencia del método.
estimacion_estratificado_proporcional <- mean(unlist(valores))
varianza_estratificado_proporcional <- valores |>
purrr::map_dbl(var) |>
weighted.mean(w = estratos$probabilidad) |>
magrittr::divide_by(n)
eficiencia_estratificado_proporcional <-
1 / (varianza_estratificado_proporcional *
coste_estratificado_proporcional)
Consideramos los mismos estratos que antes y, por tanto, la misma forma de generar valores en cada uno de ellos.
Ahora replicamos el proceso de generar en cada estrato una cantidad
óptima de valores, determinada mediante un procedimiento en dos etapas,
y aplicarles la función g a cada uno de ellos. La
estimación del coste en tiempo debe, por tanto, abarcar las dos
etapas.
n_tanteo <- 50 * cantidad_estratos
coste_estratificado_optimo <- bench::mark(
{
# Estimación de las varianzas de los estratos
n_estratos <- (n_tanteo * estratos$probabilidad) |>
ceiling() |>
pmax(2)
valores <- purrr::map2(
seq_len(cantidad_estratos),
n_estratos,
\(numero_estrato, n_estrato) {
replicate(n_estrato, {
u <- genera_vector_aleatorio_en_estrato(numero_estrato)
g(u)
})
}
)
# Cantidad óptima de valores en cada estrato
sigmas <- purrr::map_dbl(valores, sd)
n_estratos <-
(n * estratos$probabilidad * sigmas /
sum(estratos$probabilidad * sigmas)) |>
ceiling() |>
pmax(2)
# Generación de valores en cada estrato
valores <- purrr::map2(
seq_len(cantidad_estratos),
n_estratos,
\(numero_estrato, n_estrato) {
replicate(n_estrato, {
u <- genera_vector_aleatorio_en_estrato(numero_estrato)
g(u)
})
}
)
},
iterations = 10,
time_unit = unidad_de_tiempo
)$median
Finalmente, estimamos el valor de la integral, la varianza de esa estimación y la eficiencia del método.
estimacion_estratificado_optimo <- valores |>
purrr::map_dbl(mean) |>
weighted.mean(w = estratos$probabilidad)
varianza_estratificado_optimo <- valores |>
purrr::map_dbl(sd) |>
weighted.mean(w = estratos$probabilidad) |>
magrittr::raise_to_power(2) |>
magrittr::divide_by(n)
eficiencia_estratificado_optimo <-
1 / (varianza_estratificado_optimo *
coste_estratificado_optimo)
Consideramos como densidad instrumental el producto de la densidad de una \(\mathrm{Beta}(3, 1)\) para \(u_{1}\) por la densidad de una \(\mathrm{Beta}(2, 1)\) para \(u_{2}\) (tratando de esta forma a las dos coordenadas de manera independiente).
Se puede comprobar gráficamente que la densidad instrumental es «parecida» a \(g(u_{1}, u_{2}) f(u_{1}, u_{2})\).
alfa_instrumental_u_1 <- 3
beta_instrumental_u_1 <- 1
alfa_instrumental_u_2 <- 2
beta_instrumental_u_2 <- 1
densidad_nominal <- function(u) {
u_1 <- u[[1]]
u_2 <- u[[2]]
dunif(u_1) * dunif(u_2)
}
densidad_instrumental <- function(u) {
u_1 <- u[[1]]
u_2 <- u[[2]]
dbeta(u_1,
shape1 = alfa_instrumental_u_1,
shape2 = beta_instrumental_u_1) *
dbeta(u_2,
shape1 = alfa_instrumental_u_2,
shape2 = beta_instrumental_u_2)
}
library(plotly)
malla <- data.frame(
u_1 = seq_len(100) / 100,
u_2 = seq_len(100) / 100
)
malla |>
plot_ly(x = ~u_1, y = ~u_2) |>
add_surface(z = outer(
malla$u_1, malla$u_2,
Vectorize(function(u_1, u_2) {
u <- c(u_1, u_2)
g(u) * densidad_nominal(u)
})
)) |>
add_surface(z = outer(
malla$u_1, malla$u_2,
Vectorize(function(u_1, u_2) {
u <- c(u_1, u_2)
densidad_instrumental(u)
})
))
Sin embargo, en la cola de la densidad instrumental la razón de verosimilitud tiende a infinito, sin que esto quede compensado por valores pequeños de \(g\). Esto producirá una varianza infinita al aplicar el método de muestreo por importancia.
malla |>
plot_ly(x = ~u_1, y = ~u_2) |>
add_surface(z = outer(
malla$u_1, malla$u_2,
Vectorize(function(u_1, u_2) {
u <- c(u_1, u_2)
g(u)
})
)) |>
add_surface(z = outer(
malla$u_1, malla$u_2,
Vectorize(function(u_1, u_2) {
u <- c(u_1, u_2)
densidad_nominal(u) / densidad_instrumental(u)
})
))
En primer lugar, establecemos la forma de generar valores aleatorios, que ahora será a partir de la distribución instrumental escogida. Por otra parte, el producto de \(g\) por la razón de verosimilitud (es decir, el cociente entre la función de densidad nominal y la función de densidad instrumental) es el siguiente: \[\begin{align*} g(u_1, u_2) \frac{1}{\frac{u_{1}^{2}}{B(3, 1)} \frac{u_{2}}{B(2, 1)}} &= B(3, 1) B(2, 1) \frac{5 (5 u_{1} + 2)^{2} (3 u_{2} + 1)}{u_{1}^{2} u_{2}} \\ &= \frac{\Gamma(3) \Gamma(1)}{\Gamma(4)} \frac{\Gamma(2) \Gamma(1)}{\Gamma(3)} 5 \Biggl( \frac{5 u_{1} + 2}{u_{1}} \Biggr)^{2} \frac{3 u_{2} + 1}{u_2} \\ &= \frac{5}{6} \Biggl( 5 + \frac{2}{u_{1}} \Biggr)^{2} \Biggl( 3 + \frac{1}{u_2} \Biggr) \end{align*}\]
genera_vector_aleatorio <- function() {
c(
rbeta(1,
shape1 = alfa_instrumental_u_1,
shape2 = beta_instrumental_u_1),
rbeta(1,
shape1 = alfa_instrumental_u_2,
shape2 = beta_instrumental_u_2)
)
}
razon_de_verosimilitud <- function(u) {
densidad_nominal(u) / densidad_instrumental(u)
}
g_por_verosimilitud <- function(u) {
g(u) * razon_de_verosimilitud(u)
}
# Haciendo uso de la simplificación matemática indicada antes, una
# implementación más eficiente para α_u_1=3, β_u_1=1, α_u_2=2 y β_u_2=1 sería
# g_por_verosimilitud <- function(u) {
# u_1 <- u[[1]]
# u_2 <- u[[2]]
# (5 / 6) * (5 + 2 / u_1)^2 * (3 + 1 / u_2)
# }
Ahora replicamos el proceso de generar valores según la densidad
instrumental y aplicarles la función g_por_verosimilitud a
cada uno de ellos.
coste_importancia_1 <- bench::mark(
{
valores <- replicate(n, {
u <- genera_vector_aleatorio()
g_por_verosimilitud(u)
})
},
iterations = 10,
time_unit = unidad_de_tiempo
)$median
Finalmente, estimamos el valor de la integral, la varianza de esa estimación y la eficiencia del método.
estimacion_importancia_1 <- mean(valores)
varianza_importancia_1 <- var(valores) / n
eficiencia_importancia_1 <-
1 / (varianza_importancia_1 * coste_importancia_1)
Una técnica para conseguir una densidad instrumental con razón de verosimilitud acotada es perturbar ligeramente la densidad nominal. Por ejemplo, puesto que la distribución uniforme se corresponde con una distribución \(\mathrm{Beta}(1, 1)\), podemos considerar como densidad instrumental el producto de dos \(\mathrm{Beta}(1, 1/2)\).
Se puede comprobar gráficamente que la razón de verosimilitud ahora está acotada.
alfa_instrumental_u_1 <- 1
beta_instrumental_u_1 <- 1 / 2
alfa_instrumental_u_2 <- 1
beta_instrumental_u_2 <- 1 / 2
densidad_instrumental <- function(u) {
u_1 <- u[[1]]
u_2 <- u[[2]]
dbeta(u_1,
shape1 = alfa_instrumental_u_1,
shape2 = beta_instrumental_u_1) *
dbeta(u_2,
shape1 = alfa_instrumental_u_2,
shape2 = beta_instrumental_u_2)
}
malla |>
plot_ly(x = ~u_1, y = ~u_2) |>
add_surface(z = outer(
malla$u_1, malla$u_2,
Vectorize(function(u_1, u_2) {
u <- c(u_1, u_2)
densidad_nominal(u) / densidad_instrumental(u)
})
))
El producto de \(g\) por la razón de verosimilitud es ahora el siguiente: \[\begin{align*} g(u_1, u_2) \frac{1}{\frac{(1 - u_{1})^{-1/2}}{B(1, 1/2)} \frac{(1 - u_{2})^{-1/2}}{B(1, 1/2)}} &= B(1, 1/2) B(1, 1/2) \frac{5 (5 u_{1} + 2)^{2} (3 u_{2} + 1)}{(1 - u_{1})^{-1/2} (1 - u_{2})^{-1/2}} \\ &= 20 (5 u_{1} + 2)^{2} (3 u_{2} + 1) \sqrt{(1 - u_{1}) (1 - u_{2})} \end{align*}\]
genera_vector_aleatorio <- function() {
c(
rbeta(1,
shape1 = alfa_instrumental_u_1,
shape2 = beta_instrumental_u_1),
rbeta(1,
shape1 = alfa_instrumental_u_2,
shape2 = beta_instrumental_u_2)
)
}
razon_de_verosimilitud <- function(u) {
densidad_nominal(u) / densidad_instrumental(u)
}
g_por_verosimilitud <- function(u) {
g(u) * razon_de_verosimilitud(u)
}
# Haciendo uso de la simplificación matemática indicada antes, una
# implementación más eficiente para α_u_1=1, β_u_1=1/2, α_u_2=1 y β_u_2=1/2 sería
# g_por_verosimilitud <- function(u) {
# u_1 <- u[[1]]
# u_2 <- u[[2]]
# 20 * (5 * u_1 + 2)^2 * (3 * u_2 + 1) * sqrt((1 - u_1) * (1 - u_2))
# }
Ahora replicamos el proceso de generar valores según la densidad
instrumental y aplicarles la función g_por_verosimilitud a
cada uno de ellos.
coste_importancia_2 <- bench::mark(
{
valores <- replicate(n, {
u <- genera_vector_aleatorio()
g_por_verosimilitud(u)
})
},
iterations = 10,
time_unit = unidad_de_tiempo
)$median
Finalmente, estimamos el valor de la integral, la varianza de esa estimación y la eficiencia del método.
estimacion_importancia_2 <- mean(valores)
varianza_importancia_2 <- var(valores) / n
eficiencia_importancia_2 <-
1 / (varianza_importancia_2 * coste_importancia_2)
Otra técnica, conocida como muestreo por importancia defensivo, consiste en considerar como densidad instrumental una mixtura de la densidad deseada con la densidad nominal. De esta manera tenemos asegurado que la razón de verosimilitud está acotada, ya que \[\begin{equation*} \frac{f_1(x)}{p f_1(x) + (1 - p) f_2(x)} \leq \frac{f_1(x)}{p f_1(x)} \leq \frac{1}{p} \end{equation*}\]
alfa_instrumental_u_1 <- 3
beta_instrumental_u_1 <- 1
alfa_instrumental_u_2 <- 2
beta_instrumental_u_2 <- 1
densidad_instrumental <- function(u) {
u_1 <- u[[1]]
u_2 <- u[[2]]
f_1 <- densidad_nominal(u)
f_2 <-
dbeta(u_1,
shape1 = alfa_instrumental_u_1,
shape2 = beta_instrumental_u_1) *
dbeta(u_2,
shape1 = alfa_instrumental_u_2,
shape2 = beta_instrumental_u_2)
0.7 * f_1 + 0.3 * f_2
}
malla |>
plot_ly(x = ~u_1, y = ~u_2) |>
add_surface(z = outer(
malla$u_1, malla$u_2,
Vectorize(function(u_1, u_2) {
u <- c(u_1, u_2)
densidad_nominal(u) / densidad_instrumental(u)
})
))
En este caso no es posible obtener una expresión simplificada del producto de \(g\) por la razón de verosimilitud.
genera_vector_aleatorio <- function() {
if (runif(1) < 0.7) {
runif(2)
} else {
c(
rbeta(1,
shape1 = alfa_instrumental_u_1,
shape2 = beta_instrumental_u_1),
rbeta(1,
shape1 = alfa_instrumental_u_2,
shape2 = beta_instrumental_u_2)
)
}
}
razon_de_verosimilitud <- function(u) {
densidad_nominal(u) / densidad_instrumental(u)
}
g_por_verosimilitud <- function(u) {
g(u) * razon_de_verosimilitud(u)
}
Ahora replicamos el proceso de generar valores según la densidad
instrumental y aplicarles la función g_por_verosimilitud a
cada uno de ellos.
coste_importancia_3 <- bench::mark(
{
valores <- replicate(n, {
u <- genera_vector_aleatorio()
g_por_verosimilitud(u)
})
},
iterations = 10,
time_unit = unidad_de_tiempo
)$median
Finalmente, estimamos el valor de la integral, la varianza de esa estimación y la eficiencia del método.
estimacion_importancia_3 <- mean(valores)
varianza_importancia_3 <- var(valores) / n
eficiencia_importancia_3 <-
1 / (varianza_importancia_3 * coste_importancia_3)
knitr::kable(
data.frame(
`Método` = c(
"Directo",
"Estratificado proporcional",
"Estratificado óptimo",
"Importancia: versión 1",
"Importancia: versión 2",
"Importancia: versión 3"
),
`Estimación` = c(
estimacion_directo,
estimacion_estratificado_proporcional,
estimacion_estratificado_optimo,
estimacion_importancia_1,
estimacion_importancia_2,
estimacion_importancia_3
),
Varianza = c(
varianza_directo,
varianza_estratificado_proporcional,
varianza_estratificado_optimo,
varianza_importancia_1,
varianza_importancia_2,
varianza_importancia_3
),
Coste = c(
coste_directo,
coste_estratificado_proporcional,
coste_estratificado_optimo,
coste_importancia_1,
coste_importancia_2,
coste_importancia_3
),
Eficiencia = c(
eficiencia_directo,
eficiencia_estratificado_proporcional,
eficiencia_estratificado_optimo,
eficiencia_importancia_1,
eficiencia_importancia_2,
eficiencia_importancia_3
)
),
digits = 10
)
| Método | Estimación | Varianza | Coste | Eficiencia |
|---|---|---|---|---|
| Directo | 279.6914 | 3.9836253 | 0.03419515 | 7.3410305 |
| Estratificado proporcional | 279.4703 | 0.1762562 | 0.88960109 | 6.3776451 |
| Estratificado óptimo | 279.5725 | 0.1485181 | 1.03732708 | 6.4909000 |
| Importancia: versión 1 | 283.0876 | 19.0896688 | 0.15074517 | 0.3475027 |
| Importancia: versión 2 | 278.9978 | 2.8236003 | 0.17674469 | 2.0037818 |
| Importancia: versión 3 | 277.8903 | 1.0186354 | 0.20140209 | 4.8743562 |